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Abstract 

The linear strong-property-fluctuation theory (SPFT) was developed in order to estimate the constitutive 
parameters of certain homogenized composite materials (HCMs) in the long-wavelength regime. The compo- 
nent materials of the HCM were generally orthorhombic mm2 piezoelectric materials, which were randomly 
distributed as oriented ellipsoidal particles. At the second-order level of approximation, wherein a two-point 
correlation function and its associated correlation length characterize the component material distributions, 
the SPFT estimates of the HCM constitutive parameters were expressed in terms of numerically-tractable 
two-dimensional integrals. Representative numerical calculations revealed that: (i) the lowest-order SPFT 
estimates are qualitatively similar to those provided by the corresponding Mori-Tanaka homogenization 
formalism, but differences between the two estimates become more pronounced as the component particles 
become more eccentric in shape; and (ii) the second-order SPFT estimate provides a significant correction 
to the lowest-order estimate, which reflects dissipative losses due to scattering. 

Keywords: Homogenization, strong-property-fluctuation theory, metamaterials, Mori-Tanaka formalism, 
orthorhombic piezoelectric. 

1 Introduction 

Since piezoelectric materials can convert electrical energy to mechanical energy, and vice versa, they are 
of considerable technological importance. However, bulk piezoelectric materials commonly exhibit physical 
properties which render them unsuitable for particular applications. For example, certain ceramics exhibit 
strong piezoelectric properties but their weight, malleability and acoustic impedance are not suitable for 
transducer applications [1]. Accordingly, composite piezoelectric materials are often more technologically 
attractive [2]; and these can be found in a host of applications such as in transducers, sensors, actuators 
and energy harvesting devices, for example [3, 4]. Furthermore, the recent proliferation of multifunctional 
metamaterials [5] — which often take the form of homogenized composite materials (HCMs), exhibiting 
exotic constitutive properties [6] — presents interesting possibilities for piezoelectric HCMs. 

While the estimation of elastodynamic or electromagnetic constitutive parameters of HCMs is a challeng- 
ing task, especially for anisotropic HCMs, the estimation of constitutive parameters of piezoelectric HCMs 
is more challenging due to the coupling of elastodynamic and electromagnetic fields. Numerous homogeniza- 
tion formalisms have been proposed for piezoelectric HCMs, many of which build upon Eshclby's landmark 
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description of the elastodynamic response of a single ellipsoidal particle immersed in an infinite homoge- 
neous medium [7, 8]. For example, the Mori-Tanaka [9, 10, 11], self-consistent and differential approaches 
[12] — and combinations of these [13] — feature prominently in the literature. In the following we present 
a fundamentally different approach to estimating the constitutive properties of piezoelectric HCMs. based 
on the strong-property-fluctuation theory (SPFT) . A key feature of the SPFT homogenization approach — 
which distinguishes it from other more conventional approaches — is the accommodation of higher-order 
characterizations of the distributional statistics of the HCM's component materials. 

The origins of the SPFT lie in wave propagation studies for continuously random mediums [14]. It was 
later adapted to estimate the electromagnetic [15, 16, 17], acoustic [18] and elastodynamic [19] constitu- 
tive parameters of HCMs. Within the SPFT, the estimation the HCMs constitutive parameters arises by 
successive refinements to the constitutive parameters of a homogeneous comparison medium. Iterations are 
expressed in terms of correlation functions describing the spatial distributions of the component materials. In 
principle, correlation functions of arbitrarily high order may be incorporated; but, in practice, the SPFT is 
most often implemented at the second-order level of approximation, wherein a two-point correlation function 
and its associated correlation length characterize the component material distributions. 

We establish here the linear, second order SPFT appropriate to orthorhombic mm2 piezoelectric HCMs, 
arising from component materials which are randomly distributed as oriented ellipsoidal particles. The 
theoretical development builds upon the corresponding development of the orthotropic elastodynamic SPFT 
[19, 20]. A representative numerical example is used to illustrate the theory, and results are compared with 
those from the well-established Mori-Tanaka formalism. 



2 Theory 

2.1 Preliminaries 

In the following, we consider piezoelectric materials described by constitutive relations of the form [21, 22] 



Da — (iarnnSrnn ^~ ^an-^n 



(1) 



wherein the elastic strain Smn and the electric field En are taken as independent variables, which are related 
to the stress aab and dielectric displacement Da via the elastic stiffness tensor Cabmn (measured in a constant 
electric field), the piezoelectric tensor Cnab (measured at a constant strain or electric field), and the dielectric 
tensor Can (measured at a constant strain). Here, and hereafter, tensors are represented in plain font and 
lowercase tensor indexes range from 1 to 3, with a repeated index implying summation. 

We develop the SPFT in the frequency domain. Accordingly the complex valued representations of the 
stress, strain and electromagnetic fields have an implicit exp {—iujt) dependency on time t, with uj being the 
angular frequency and i = The possibility of dissipative behaviour is thereby accommodated via the 

imaginary parts of complex- valued constitutive parameters. 

The constitutive relations (1) are more conveniently expressed in the symbolic form 



^aB — ^aBMn^Mn^ 



where the extended stress symbol 



the extended stiffness symbol 
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lb': 
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B = 6= 1,2,3 
B = 4 



(2) 
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B = 6= 1,2,3; Af = m = 1,2.3 
B = 6 = 1,2,3; M = 4 
B = 4; M = m = 1,2,3 
B,M = 4 



(4) 



2 



and the extended strain symbol 



m 



■'Mn 



= 1,2,3 
M = 4 



(5) 



Here, and hereafter, uppercase indexes range from 1 to 4. Note that the extended quantities defined in 
eqs. (3), (4) and (5) are not tensors — these are simply symbols which are introduced to allow a compact 
representation of the piezoelectric constitutive relations [10]. 

In developing the SPFT appropriate to piezoelectric HCMs, it is expedient to express the constitutive 
relations (2) in matrix- vector form as 

a=QS, (6) 

wherein a and S are column 12-vectors representing the extended stress and extended strain symbols, re- 
spectively, and C is a 12x12 matrix which represents the extended stiffness symbol. Here, and hereafter, 

matrixes arc denoted by double underlining and bold font, whereas vectors are in bold font with no un- 
derlining. For use later on, we note that the pqth entry of a matrix A is written as [ A]^^, while the pth 

entry of a vector b is written as [h]^. Accordingly, the matrix entry [ A • B]^^ = [ A]^^ [B ]^^, the vector 

entry [A-b]p= [ A]^^ [b] ^, and the scalar a • b = [aj^ [bj^. The adjoint, determinant, inverse, trace and 

transpose of a matrix A are denoted by adj ( A), det ( A), A~^, tr ( A) and A'^, respectively. The n x n 
null matrix is written as 

= nxn 

Our concern in this article is with orthorhombic mm2 piezoelectric materials [21, 22]. For this symmetry 
class, the extended stiffness matrix has the block matrix form 




(7) 



where the 9x9 stiffness matrix C may be expressed as 

/ C 

= ( 

0. 



= 3x3 
V 23X3 



with the 3x3 symmetric matrix components 
C 



Cii Ci2 Ci3 
= I Ci2 C22 C23 
Ci3 C23 C33 




while the 9x3 piezoelectric matrix e may be expressed as 



(8) 
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(10) 



and the 3x3 dielectric matrix e as 



£11 
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633 



(11) 



The correspondence between the tensor /extended symbol representation and the matrix-vector representa- 
tion is described in Appendix A. 
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In an analogous fashion, the material density p may be represented via the extended density symbol 

r p, B = M=h2,3 
" I 0, otherwise ' ^^^> 

which has the 4x4 extended matrix counterpart p with entries 



P 

= Jmp 



Pmp- (13) 



2.2 Component materials 



We consider the homogenization of a composite comprising two component materials, labelled as component 
material '1' and component material '2'. In general, both components are homogeneous, orthorhombic mm2 
piezoelectric materials, characterized by the stiffness tensors C^^mn' piezoelectric tensors e^^jj,, dielectric 
tensors ein and densities p^^^ (£ = 1,2). In conformity with the notational practices introduced in §2.1. the 
component materials are also described by the extended stiffness symbols C^^j^^^ (and their 12x12 matrix 

equivalents C*' ^) and extended density symbols p^M (^'^'^ their 4x4 matrix equivalents p*-^^). 

The component materials are randomly distributed as identically-oriented, conformal, ellipsoidal parti- 
cles. The principal axes of the ellipsoidal particles arc aligned with the Cartesian axes. Thus, the surface of 
each ellipsoidal particle may be parameterized by the vector 

r(")=r/n-f, (14) 

where ry is a linear measure of size, f is the radial unit vector and the diagonal shape matrix 



U = ^^ 6 0, (a,6,ceM+). (15) 




Let V denote the space occupied by the composite material. Then V = V^^^ U V^"^^ where V^^^ and 

V^'^'' contain the two component materials labelled as '1' and '2', respectively, and V^^^ n V^"^^ = 0. The 
distributional statistics of the component materials are described in terms of moments of the characteristic 
functions 

r 1, reFW, 

$(^)(r) = < (^ = 1,2). (16) 

I 0, r^VW, 

The first statistical moment of i.e., 

($W(r))=/(^), (^=1,2), (17) 

2 

delivers the volume fraction of component material £, which is subject to the constraint ^^/^^^ = 1- The 

second statistical moment of constitutes a two point covariance function. Investigations involving 
the electromagnetic SPFT have demonstrated that the specific form of the covariance function has only a 
minor influence on estimates of HCM constitutive parameters, for a range of physically-plausible covariance 
functions [25] . Here we adopt the physically-motivated form [26] 

r ( $(^) (r) ) ($ W (r') ) , I • (r - r') I > i 

($W(r)$W(r')) = < , (18) 

I (<I>('Hr)), \g-'-ir-r')\<L 

which has been widely used in electromagnetic and elastodynamic SPFT studies. The correlation length L 
in eq. (18) is required to be much smaller than the associated piezoelectric wavelengths, but larger than the 
particle size parameter r/. 
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2.3 Comparison material 

A homogeneous comparison material provides the initial ansatz for an iterative procedure that delivers a 
succession of SPFT estimates of the HCM constitutive parameters [19] . Accordingly, the comparison material 
represents the lowest-order SPFT estimate of the HCM. In consonance with the component materials, the 

comparison material is an orthorhombic mm2 piezoelectric material, in general. The piezoelectric constitutive 
properties of this orthorhombic comparison material (OCM) arc encapsulated by its extended stiffness symbol 

^iMi'q (^md its 12x12 matrix equivalent C ) and extended density symbol p^p"' (and its 4x4 matrix 

equivalent 

In order to establish the spectral Green function for the OCM — which is a key element in the SPFT 
formulation — we first consider the corresponding extended equation of motion. This may be written in the 
frequency domain as [27] 

^iMpldidqUp + w^um = -Fm, (19) 

where the extended displacement 

{ Um, M = m = 1,2,3 , . 

= I p = 4 ^^^> 

combines the displacement Um and electric scalar potential i>, and the extended body force 

^ ( F^, M = m = 1,2,3 , . 

= I -g, M = 4 (21) 

combines the body force F^ and the electric charge q. Accordingly, the sought after spectral Green function 
for the OCM emerges as the 4x4 matrix 

g(ocm) (1^) ^ [fc2^(k) - ]-\ (22) 

where the 4x4 matrix a(k) has entries 

a(k) = (23) 

- J MP fc^ 

Herein, k = fck = (fci, ^2, fcs) with k = {sin 9 cos cj), sin 61 sin cos 6). 

The OCM extended constitutive symbols C^^p^ and P^mp^^ ^'^^ derived via the imposition of the two 
conditions [19, eqs. (2.72), (2.73)] 

($W(r)^Wp^ + $(2)(r)^l^p,) = 0, (24) 
($(^) (r) - p^"^")] + $(2) (r) [^(2) _ p{ocm)-\ ^ ^ 0^ ^25) 



MP 



MP 



which is necessary to remove certain secular terms. In eq. (24), the quantities 

^[mPq — ipiMSt ~ ^IMSt^ ''l^StPq-' ~ ^,'^)-, (26) 



where viglp^ is given implicitly through 



'^Pq — nPqStJSt ' l^' J 

iTj - ^Tj + ^TjlM y^lMPq " '^IMPq) ^Pq' l^^) 
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with the renormalization tensor 



d9 



smf 



— X 



WpstU = < 



8^ Jo '^'^ Jo (H"' • k) • (H"' • k) 

(n-^-k)t{(n-'-k)«[i-'(n-^-k)],^ 

+(n-'.k),[^-i(n-^.k)]^^}, 



27T pTT 



— d4) dOsmO 
Stt Jo Jo 



(n-'-k)t(n-'-k)«[i-'(n-'-k)] 



p(7 



(u-^.k).(u-^.k) 



P = p = 1,2,3 • 



P =4 



(29) 



Upon substituting eqs. (26)-(28) into eq. (24), exploiting eq. (17), and after some algebraic manipula- 
tions, we obtain 



/ 



(1) 



/ 



(2) 



: 12x12' 



(30) 



wherein the 12x12 matrix equivalent of Wustu (namely, W) has been introduced and ^ denotes the matrix 
operation defined in Appendix A. The OCM stiffness matrix may be extracted from (30) as 



C^"+/(') T + (C^"'-C™)-W 



(2) iL(ocm). 



(31) 



where r is the 12 x 12 matrix representation of the extended identity TrSTu = TRstu, as described in 
Appendix A. By standard numerical procedures, such as the Jacobi method [28], the nonlinear relation (31) 

is solved for c^"*^™' . 

After combining eq. (17) with cq. (25), it follows immediately that the OCM density is the volume 
average of the densities of the component materials '1' and '2'; i.e.. 



^(ocm) 



/(l)^(l) +/(2)^(2)_ 



(32) 



2.4 Second-order SPFT 

Building upon the corresponding results for the elastodynamic SPFT [19], the second-order^ estimates of 
the HCM extended stiffness and density symbols may be expressed in terms of three-dimensional integrals 
as 



XY 



G(°'=")(k) 



YU 



{ks [i"'(k) 



(33) 



XY 



and 



PmP^ = PmT^ +^'J d'k BMSUp{k) (k) 



su 



(34) 



"^The first-order SPFT estimate is identical to the zeroth-order SPFT estimate which is represented by the comparison 
material. 
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The symbols B\.^p^{]s) and Bmsup{}^ represent the spectral covariance functions given as 

(2) _ ^(1) ^ _ Ai 



^tUPq W - 



1) 
nUPq 



j dPR r(R) exp {-ik ■ R) 



with 



, , \Pms~ Pms) \Pup~ Pup) f ^ , , 
Bmsup{^ = ^ j d^RT{K) exp(-ik. R) 



r(R)=r(r-r') = ($(i)(r)$(i)(r')) - ($^^^(r)) ($(i)(r')) 
= ($(2)(r)$(2)(r')) _ ($(2)(r)) ($(2)(r')). 



(35) 



(36) 

In order to make the integrals in the expressions for cj^pg and p^^p*^ presented in eqs. (33) and (34) 
numerically tractable, we simplify them as follows. Let us begin with the integral on the right sides of eqs. 
(35). Upon implementing the step function-shaped covariance function (18), wc find 



/ d^R r(R) exp {-ik • R) = / d^R exp [-i (H • k) • R] . 

J •'|R|<i' 

Thereby, the expressions for B\ffp^ (k) and Bmsup (k) reduce to 
Hi)fi2)(A2) _Ai) \(A2) _Ai) \ 

J J yilMRs ^IMRs J yitUPq ^tUPq J sm [kaL) 



(37) 



r/lMRs n.\ 



BmsupO^) 



' ■ ■ - 



2 {TTkaY 



ka 
sin {kaL) 



L cos {kaL) 



ka 



— L cos {kaL) 



(38) 



wherein the scalar function 

a = a{6, 4>) = ^ a? s\v? 6 cos^ (j) + h'^ s\v? 9 s\v? (j) + c^ cos^ 6 



(39) 



is introduced. 

Now we turn to the integrals in (33) and (34). In analogy with the corresponding expression for the 
elastodynamic SPFT [20] , the spectral Green function G^"'^'"^ (k) may be conveniently expressed as 



where the 4x4 matrix function 
and the scalar function 



I ^ D(k) 
= ~ A(k) ' 

D(k) = adj [ A;2 a(k) - w^p^"^™) 



(40) 
(41) 



A{k) = k^det [^(k)] - tr {adj [fc'i(k)] • co^p^"^"''^} - fc^tr [adj(a;2p('"=")) • ^(k) 



•'('■■{ 


a 


J 4' 


^ [^«(k)-adj(a 






a(k) adj {uy 


l(£) 


42 


i(k) 


24 


adj(o.V) 


22 


i(k) 


31 


i(k) 


13 


adj(u;V) 



11 



33 . 



(42) 



with the 3x3 matrixes a" and p" having entries 



pq 



{p,q= 1,2,3). 



(43) 
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Through exploiting oqs. (38) and (40). the integrals in cqs. (33) and (34) with respect to k can be 
evaluated by means of calculus of residues: The roots of A(k) = give rise to seven poles in the complex-A: 
plane, located at A; = 0, ±pi, ±P2, ips, which are chosen such that p„ (n = 1,2,3) lie in the upper half of 
the complex plane. Prom eq. (42), we find that the nonzero poles satisfy 



pI 



1 /2V3PB Pc 



Pa 



Pa + 



1 / (l + iV3)PB _ {l-iV2.)Pc \ 
3 1^ 22/3Pc 24/3 j 

1 / (1 - i^)PB _ {l+i^)Pc \ 



22/3p^ 



24/3 ) 



wherein 



Pa 

Pb 
Pc 
Pd 



w^tr |adj 


a 


(k) 




p(ocm) 1 


3 det 







—G\ + 3 C_B, 



1/3 



Pu + {^Pl + Plf\ 

-2C\ + 9CaCb-27Cc 



with 



Ca = 



—fjp' tr |adj 


i(k)^ 




det 

_ / 


i(k)^ 
i(k) 


tr [a 

J 44 L — 


det 







(44) 
(45) 

(46) 

(47) 

(48) 
(49) 
(50) 

(51) 



+ 



(k)]^J#) 



^«(k)-adj (p")] + [i(k)]^^ [=(^^]i4 H (p^°""^) 

adj(^(— ))] +[a(k)l [a(k)l [adj (p(<«='") ) 
\= /J 22 L— J43 L— J34 L \= / 



Cc 



— tr |adj 


E 


•i(k)} 


det 


i(k) 





(52) 
(53) 



Thus, by application of the Cauchy residue theorem [29], the SPFT estimates are delivered in terms of 
two-dimensional integrals as 



^IMPq 



^IMPq "I" 



47ri 



TT /•IT 



kf sin 6 



=0 J 0=0 



(ka) det a(k) 



[I' 



(ocm) 



XY 



b(k) 



YU 



and 



Pmp ~ Pmp 



+ 



{^lm4s ~ ^lm4s} {^s [l ^ (k) } j (dc/Pg ~ ^tUPq) 



2m 



2-7T p7T 

/ d(pd9 
,=0 J 0=0 det 



smf 



a(k) 



(54) 



Mk)J^^, (55) 
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where the 4x4 matrix 



b(k) 



1 [ e'^- 
2i I apf{p 



e'L'yp3-D{p3\l ■ k) 



1 — iLapi 



1 

2 2 2 

(^PiPipi 



pI){pI 



Pi) 



1 — iLaps 



+ 



1 



Pi P2 



Pi 



0-P2(Pl-P2)(P2-P3) 



1 5' ... ^■ 



1 — iLap2 



(56) 



The expressions for the second order SPFT estimates Ci^pq ^^'^ P^mp^' ^1^- i^^) ^^'^ i^^) be 
evaluated by standard numerical methods [30] . 

It is particularly noteworthy that CiMPq ^^'^ Pai'p^'' complex-valued for L > 0, even when the 

corresponding quantities for the component materials, i.e., C/MPq ^^'^ Pmp — 1>2), are real-valued. This 
reflects the fact that the SPFT accommodates losses due to scattering [17]. From energy considerations, the 
imaginary part of the extended compliance matrix, namely [22] 



Uspft) 



M 



(spft)) 



I 



\ 



{spft) 



Aspft) 



))^ • (^c'^'Pf'^ 



\ 



.(spft) 



(57) 



is required to be positive definite for passive materials [31]. The constitutive matrixes C^^^^*\ ei^pft) ^^^^ 
^{spft) ^.-^^ right side of eq. (57) are related to the extended stiffness matrix ^'^ (and thereby to the 
extended stiffness symbol c['^p^) per eq. (7). 



3 Numerical results 
3.1 Preliminaries 

In order to illustrate the theory presented in §2, let us now consider a representative numerical example. A 
comparison for the SPFT estimate of the HCM constitutive parameters is provided by the corresponding 
results computed using the Mori-Tanaka formalism [9, 12, 23, 24]. In the case of orthorhombic mm2 
piezoelectric component materials, the Mori-Tanaka estimate of the extended stiffness matrix for the HCM 
is given by [13] 



where the 12x12 matrix 



g(MT) ^ 



^{Esh) 



(MT) 



(£'")■( 



.(2) 



.(1) 



(58) 



(59) 



with S'-^*'*^ being the 12 x 12 matrix representation of the Eshelby tensor [7, 10, 32]. Details on evaluating 
^^^"''^ can be found in Appendix B. 

In the following, we present the numerical evaluation of the 12x 12 extended stiffness matrix of the HCM, 

\j {^fx CT7X ^ 

namely C , as estimated by the lowest-order SPFT (i.e., hem = oem), the second-order SPFT (i.e., 
hem = spft) and the Mori-Tanaka formalism (i.e., hem = MT). The matrix c^^""^^ has the form represented 
in eq. (7). The second-order SPFT density tensor pf^p^ ^^^'^ evaluated; the numerical evaluation of the 
lowest-order SPFT density p'mp^^ need not be presented here as this quantity is simply the volume average 



9 



of the densities of the component materials. An angular frequency of a; = 27r x 10^ was selected for all 
second-order SPFT computations. 

The eccentricities of the ellipsoidal component particles are specified by the shape parameters {a,b,c}, per 
eqs. (14) and (15). To allow direct comparison with results from previous studies [13], component material 
'1' was taken to be the piezoelectric material polyvinylidene fluoride (PVDF) while component material '2' 
was taken to be the thermoplastic polyimide LaRC-SI, which has no piezoelectric properties. The stiff'ness 
constitutive parameters of the component materials are tabulated in Table 1. The nonzero piezoelectric 
constitutive parameters of PVDF are: ens = esi = 0.024, 6223 = 632 = 0.001 and 6333 = 633 = —0.027 in 
units of C m~^. The dielectric constitutive parameters of PVDF are: en = 7.4, £22 = 9.6 and £33 = 7.6, 
whereas those of LaRC-Sl arc: en = £22 = £33 = 2.8, all in units of e„ = 8.854 x 10^"'^^ F m"-'^ (the 
permittivity of free space). Lastly, the densities of PVDF and LaRC-SI are 1750 and 1376, respectively, in 
units of kg m~^. 



Stiffness parameter 


PVDF (GPa) 


LaRC-SI (GPa) 


Ciiii = C'li 


3.8 


8.1 


^1122 — ^12 


1.9 


5.4 


C1133 = Ci3 


1.0 


5.4 


^2222 = ^22 


3.2 


8.1 


^2233 = ^23 


0.9 


5.4 


C3333 = C'33 


1.2 


8.1 


^2323 = ^44 


0.7 


1.4 


Cl313 = C'55 


0.9 


1.4 


C1212 = C'ee 


0.9 


1.4 



Table 1: The stifi^ness constitutive parameters of the component materials in units of GPa (after [13]). 



3.2 Lowest-order SPFT 

We begin by considering the lowest-order SPFT estimates of the HCM constitutive parameters. In Fig. 1, 
components of the HCM extended stiffness matrix C , as computed using the lowest-order SPFT and the 
Mori-Tanaka formalism, are plotted as functions of volume fraction /^^^ for the case where the component 

particles are spherical (i.e., a = b = c). Plots of only a representative selection of the components of C 
are presented in Fig. 1; plots for those components which are not presented in Fig. 1 are qualitatively similar 
to those that are presented. Only relatively minor differences between the lowest-order SPFT estimates and 
the Mori-Tanaka estimates are observed, with the differences between the two being greatest for mid-range 
values of /^^^ . Plots for both the SPFT and Mori-Tanaka estimates are necessarily constrained by the limits 

hm c''-' = &'\ lim C^'^-^=c''\ (60) 

I ' ■ corresponding graphs for the cases where the components particles are described by the shape param- 
eters {a/c — 5, b/c — 1.5} and {a/c = 10, b/c = 2} are provided in Figs. 2 and 3, respectively. A comparison 
of Figs. 1-3 reveals that the differences between the lowest-order SPFT and Mori-Tanaka estimates are 
accentuated as the component particles become more eccentric in shape, especially at mid-range values of 
/(^) for the piezoelectric parameters and the dielectric parameters. 

3.3 Second-order SPFT estimate 

Now let us turn to the second-order SPFT estimates of the HCM constitutive parameters. We considered 
these quantities as functions of kL, where k is an approximate upper bound on the wavenumbers supported 
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by the HCM, as estimated by [20] 




(61) 



wherein 



A = - 



-t 

=1 



M = 



(|[fi' 



J 12 



2 



+ 



c(^) 



C(^) 



13 



55 



+ 



c(^) 

— J 23 

c(^) 

— J 66 



(62) 



and L is the correlation length associated with the the two-point covariance function (18). In Fig. 4, the real 



~ (sp/t) 

and imaginary parts of the components of C 



C — C are plotted against kL for ' = 0.5. 



The values of the shape parameters {a, b, c} correspond to those used in the calculations for Figs. 1-3. As in 

~ (spft) 

§3.2, only a representative selection of the components of C are plotted in Fig. 4; the graphs for those 

components that are not represented in Fig. 4 are qnalitatively similar to the graphs which do appear. 

The second-order corrections to the lowest-order SPFT estimates are observed in Fig. 4 to grow expo- 
nentially in magnitude as the correlation length increases from zero. Furthermore, the magnitudes of both 

{spft) 

the real and imaginary parts of C generally grow faster with increasingly correlation length when the 
components particles are more eccentric in shape. At L = 0, the second-order and lowest-order SPFT 
estimates coincide. While the second-order corrections are relatively small compared to the lowest-order 
SPFT estimates, a highly significant feature of the second-order corrections is that these are complex-valued 

are purely real-valued. We note that for all 



. . ^("-^b) "(ocm 

With nonzero imaginary parts, even though C and G 



computations the imaginary part of the extended compliance matrix M was found to be positive def- 

^ [spft) 

inite, which corresponds to positive loss [31]. Thus, the emergence of nonzero imaginary parts of C 
indicates that the HCM has acquired a dissipative nature, despite the component materials being nondissi- 
pative. The dissipation is attributed to scattering losses, since the second-order SPFT takes into account 
interactions between spatially-distinct scattering particles via the two-point covariance function (18). As 
the correlation length increases, the number of scattering particles that can mutually interact also increases, 
thereby increasing the scattering loss per unit volume. 

Finally, we turn to the second-order SPFT estimate of the HCM density. The real and imaginary 



parts of the matrix entry 



wherein p^'^P/*) = p' 



p^""""', are plotted as functions of kL in 



and 



22 



(spft) _ p(ocm 

are much the same as those for 



33 



-Sspft) 



11 



Fig. 5. The corresponding graphs for 

but with minor differences in magnitudes. The second-order SPFT estimates of the HCM density exhibit 
characteristics similar to those of the corresponding HCM stiffness, piezoelectric and dielectric constitutive 



parameters. That is, lim^pi*/^*^ = p^'"^"^ and 



.(spft) 



Jocm) 



for a = 1,2 and 3. Also, the differences 



between pf^'^Pf^^ and p^'"^'") increase exponentially as the correlation length increases, and this effect is most 
accentuated when the component particles are most eccentric in shape. 

We remark that a complex-valued, anisotropic density also crops up in the second-order elastodynamic 
SPFT for orthotropic HCMs [20], as well as in other homogenization scenarios [33, 34]. 



4 Closing remarks 

The linear SPFT has been fully developed for the case of orthorhombic mm2 piezoelectric HCMs, based on 
component materials distributed as oriented ellipsoidal particles. The multifunctionality of such HCMs is 
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central to the notion of mctamatcrials [5] . The second order estimates of the HCM constitutive parameters 
are expressed in terms of numerically-tractable two-dimensional integrals, for a specific choice of two-point 
covariance function. This theoretical result further extends the application of the SPFT in the homogeniza- 
tion of complex composites, effectively bridging the clastodynamic SPFT for orthotropic HCMs [19, 20] and 
the electromagnetic SPFT for anisotropic dielectric HCMs [35, 36]. Furthermore, the path has now been 
cleared towards the development of the SPFT for piezoelectric/piezomagnetic HCMs [37], with bianisotropic 
electromagnetic properties [17]. Let us remark that the mathematical description of piezoelectric HCMs 
presented herein also extends to electrokinetic processes [38]. 

From our theoretical considerations and representative numerical studies, the following conclusions were 
drawn: 

• The lowest order SPFT estimate of the stiffness, piezoelectric and dielectric properties of the HCM 
are qualitatively similar to those estimates provided by the Mori-Tanaka formalism. 

• Differences between the estimates of the lowest order SPFT and the Mori-Tanaka formalism are great- 
est at mid-range values of the volume fraction, and accentuated when the component particles are 
eccentric in shape. 

• The second-order SPFT provides a correction to the lowest-order estimate of the HCM constitutive 

properties. The magnitude of this correction is generally larger when the component particles are more 
eccentric in shape, and vanishes as the correlation length tends to zero. 

• While the correction provided by the second-order SPFT is relatively small in magnitude, it is highly 
significant as it indicates dissipation due to scattering loss. 



Appendix A 



The extended symbol AaMPq {o-Tq £ {1,2,3}, M,P e {1,2,3,4}) may be conveniently represented by the 

(7, K e [1, 12]), upon replacing the index pair aM with 7 and the index 



12 X 12 matrix with entries 



pair Pq with k. For the most general 12x 12 matrix encountered in this paper, which has the form 



A = 



/ 




Al,2 


-4l,3 


























-4l,12 




A2,l 


A2,2 


-42,3 


























-42,12 




As,! 


-43,2 


-43,3 


























-43,12 













-44,4 








-44,4 











-44,11 



















-45,5 








-45,5 





-45,10 

























-46,6 








-46,6 






















-44,4 








-44,4 











-44,11 



















-45,5 








A5,5 





-45,10 

























-46,6 








-46,6 

























-4l0,5 








-4l0,5 





-4io,io 



















-4ll,4 








-4ll,4 











-4ii,ii 





\ 


-4l2,l 


-4.12,2 


-4l2,3 


























-4l2,12 



(63) 



/ 

the correspondence between the extended symbol indexes and the matrix indexes is provided in Table 2. The 
scheme presented in Table 2 also relates the extended symbol taM to the corresponding column 12-vector 
entries [t] . 

" t 

We introduce the matrix A which plays a role similar to the matrix inverse insofar as 



A^ • A = A • A^ = T. 



(64) 
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aM or Pq 


7 or K 


aM or Pq 


7 or K 


aM or Pq 


7 or K 


aM or Pq 


7 or «; 


11 


1 


23 or 32 


4 


23 or 32 


7 


14 or 41 


10 


22 


2 


13 or 31 


5 


13 or 31 


8 


24 or 42 


11 


33 


3 


12 or 21 


6 


12 or 21 


9 


34 or 43 


12 



Table 2: Conversion between extended symbol and matrix notation. 



Herein, 



/ 1 


^3x3 


^3x3 


^3x3 


= 3x3 






^3x3 


§3X3 






23X3 


V §3X3 


= 3x3 


= 3x3 


I 



(65) 



is the 12x12 matrix representation of the extended identity symbol, with I being the 3x3 identity matrix, 
and we have 

A • r = r • A = A- (66) 

" t 

The matrix A has the form 



tl,l 


tl,2 


tl,3 


























tl,12 


t2,l 


t2,2 


t2,3 


























t2,12 


t3,l 


t3,2 


t3,3 


























t3a2 











t4,4 

2 








t4,4 

2 











t4,ll 

















t5,5 

2 








t5,5 

2 





ts.io 























t6,6 

2 








te.e 
2 




















t4,4 

2 








t4,4 

2 











t4,ll 

















ts.s 
2 








t5,5 

2 





ts.io 























t6,6 

2 








2 























tlO.5 








tl0,5 





tio.io 

















til, 4 








tll,4 











til, 11 





tl2,l 


tl2,2 


tl2,3 


























tl2,12 



(67) 
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with entries 



tl,l — (—^12,3^2,2^3,12 + ^12,2^2,3^3,12 + ^12,3^2,12^3,2 — ^12,12^2,3^3,2 — 

-4l2,2-42,12^3,3 + ^12,12^2,2-43,3)/A, (68) 
tl,2 = (^1,2^12,3^3,12 — ^12,2^1.3^3,12 — ^1,12^12,3^3,2 + ^12,12^1,3^3,2 — 

^1,2^12,12^3,3 + ^1,12^12,2^3,3)/A, (69) 
tl,3 = (—^1,2^12,3^2,12 + ^12,2^1,3^2,12 + ^1,12^12,3^2,2 — ^12,12^1,3^2,2 + 

vll,2Ai2,12^2,3 - ^l,12Ai2,2^2,3)/A, (70) 
t2,l = (—^12,3^2,12^3,1 + ^12,12^2,3^3,1 + ^12,3^2,1^3,12 — ^12,1^2,3^3,12 — 

-4l2,12^2,l-43,3 + ^12,l-42,12-43,3)/A, (71) 
t2,2 = (^1,12^12,3^3,1 — ^12,12^1,3^3,1 — ^1,1^12,3^3,12 + ^12,1^1,3^3,12 — 

^1,12^12,1^3,3 + ^1,1^12,12A3,3)/A, (72) 
t2,3 = (—^1,12^12,3^2,1 + ^12,12^1,3^2,1 + ^1,1^12,3^2,12 — ^12,1^1,3^2,12 + 

^1,12^12,1^2,3 - ^I,l>ll2,12^2,3)/A, (73) 
ts.l = (^12,2^2,12^3,1 — ^12,12^2,2^3,1 — ^12,2^2,1^3,12 + ^12,1^2,2^3,12 + 

^12,12^2,1^3,2 - ^12,1^2,12^3,2)/A, (74) 
t3,2 = (^1,2^12,12^3,1 — ^1,12^12,2^3,1 — ^1,2^12,1^3,12 + ^1,1^12,2^3,12 + 

vll,12Ai2,l A3,2 - Al,lAl2,12A3,2)/A, (75) 
t3,3 = (—^1,2^12,12^2,1 + ^1,12^12,2^2,1 + ^1,2^12,1^2,12 — ^1,1^12,2^2,12 — 

-41,12^12,1^2,2 + Al,lAl2,12^2,2)/A, (76) 

^^'^ 2(^11,11^4,4 — ^4,11^11,4)' ^ ^ 

'^'5 ~ 971 A A A T' ^'^> 

^(,^10,10^5,5 — ^5, 10^10, 5 j 

t6,6 = (79) 

tlO,10 = 77 . ^^'^A A ^' 

(,^10,10^5,5 — ^10,5^5, 10 J 

^ 71 A ^ A \' ^^^^ 

(/111, 11/14,4 — /lll,4/l4,llj 

tl2,12 = (—^1,3^2,2^3,1+^1,2^2,3^3,1+^1,3^2,1^3,2 — ^1,1^2,3^3,2 — 

^I,2>l2,l>l3,3 + Ai,iA2,2A3,3)/A, (82) 
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tl,12 
t2,12 
t3,12 

t4,ll 

ts.io 

tl2,l 
tl2,2 
tl2,3 

tll,4 
tl0,5 



(^1,3^2,2^3,12 — ^1,2^2,3^3,12 — ^1,3^2,12^3,2 + ^1,12^2,3^3,2 + 

^1,2^2,12^3,3 - Ai,i2A2,2A3,3)/A. 

(^1,3^2,12^3,1 ^ ^1,12^2,3^3,1 ^ ^1,3^2,1^3,12 +^1,1^2,3^3,12 + 
^1,12^2,1^3,3 - -414^2, 12A3,3)/A, 

(—^1,2^2,12^3,1 + ^1,12^2,2^3,1 + ^1,2^2,1^3,12 " ^1,1^2,2^3,12 — 
-4l,12^2,l-43,2 + ^l,lA2,i2^3,2)/A, 

^4^11 

2(^11,4^4,11 — ^11,11^4,4) ' 

^5,10 

2(^5,10^10,5 — ^10,10^5,5) 

(^12,3^2,2^3,1 — ^12. 2-42, .3^34 — ^12,3^2,1^3,2 + ^12,1^2,3^3,2 + 
^12,2^2,1^3,3 - Al2,lA2,2^3,3)/A, 

(—Ai. 2^12,3^3,1 + ^12.2^1,3^3,1 + ^1,1^12,3^3,2 — ^12,1^1,3^3,2 + 
^1,2^12,1^3,3 - ^1,1^12,2^3,3)/A, 

(^1,2^12,3^2,1 — ^12,2^1,3^2,1 — ^1,1^12,3^2,2 +^12,1^1,3^2,2 — 
-4l,2-4i2,lA2,3 + ^l,lAi2,2A2,3)/A, 

^11,4 

2(>lll,4^4,ll — ^11,11^4,4) ' 

^10,5 

1{A 



10,5^5,10 



110,10^5,5J 



(83) 

(84) 

(85) 
(86) 

(87) 

(88) 

(89) 

(90) 
(91) 

(92) 
(93) 



where the scalar 
A 



= Ai, 12^12,3^2, 2^3,1 — ^12,12^1,3^2,2^3,1 — ^1,1^12,3^2,2^3,12 + ^12,1^1,3^2,2^3,12 — 

^1,12^12,3^2,1^3,2 + ^12,12^1,3^2,1^3,2 + ^1.1^12,3^2,12^3,2 



^12,1^1,3^2,12^3,2 + 

^1,12^12,1^2,3^3,2 — ^1,1^12,12^2,3^3,2 — ^1,12^12,1^2,2^3,3 + ^1,1^12,12^2,2^3,3 + 
^12,2(^1,3^2,12^3,1 ^ ^1,12^2,3^3,1 ^ ^1,3^2,1^3,12 + ^1,1^2,3^3,12 + 
^1,12^2,1^3,3 — ^1,1^2,12^3,3) + ^1,2(— ^12,3^2,12^3,1 + ^12,12^2,3^3,1 + 
^12,3^2,1^3,12 — ^12,1^2,3^3,12 — ^12,12^2,1^3,3 + ^12,1^2,12^3,3)- 



(94) 



Appendix B 



The extended Eshelby symbol appropriate to orthorhombic mm2 piezoelectric materials, distributed as 
ellipsoidal particles with shape parameters {a, 6, c}, is given by [10, 32] 



S 



(esh) 
MnAh 



1 /^(l) 




<.27r 

/ dw [F^j,^(d)+Fnj,^{^)], 

h 


M = m = 1, 2, 3 


< 

1 ^(1) 

. 47r^^^^ 










/ du! F4jsn{'&) , 

/o 


M = 4 



(95) 



wherein 



"a Cl ^ C2 

i?i = — , ^2 = -r, 
a 

Ci = (l-C3')'/'cos(a;), 



KjR = ^sCi%Jn 
c 

C2 = (l-C3)'/'sin(a;), C3 = C3 



(96) 
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The integrals in eqs. (95) can be evaluated using standard numerical methods [30]. 

(esh) 

The conversion from the extended Eshclby symbol S)^^^^^^ to the extended Eshelby 12x 12 matrix, namely 
g(£;s/i)^ follows the scheme described in Appendix A. 
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Figure 1: Plots of 



' (hem) 



(in GPa), 



' {hem) 



(in C ni-2) and (l/e,,) 



(hem) 



as estimated 



using the lowest-order SPFT (i.e., hem = ocm) (black, dashed curves) and the Mori-Tanaka formalism (i.e., 
hem — MT) (red, solid curves), versus the volume fraction of component material '2'. Component material 
'1' is PVDF and component material '2' is LaRC-SI, as described in §3.1. The component materials are 
distributed as spheres (i.e., a = b = c). 
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Figure 2: As Fig. 1 but with the component materials distributed as ehipsoids with (a/c = 5 and h/c = 1.5). 
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Figure 3: As Fig. 1 but with the component materials distributed as ellipsoids with (a/c = 10 and h/c = 2). 
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Figure 4: Plots of the real and imaginary parts of the second-order SPFT estimates 



(in GPa), 



1.12 



(in C m-2) and (lOVco 



where C 



(spft) 



12,12 



^{spft) ^ (ocm) - 

C — C , versus kL, with 



/(^) = 0.5. The results from the spherical particle (i.e., a = 6 = c = 1) case (red, solid line) are plotted 
alongside the cases with elliptical particles with a — 5, b — 1.5, c — 1 (blue, short-dashed line) and a — 10, 
b = 2, c = 1 (black, long-dashed line). Component material '1' is PVDF and component material '2' is 
LaRC-SI, as described in §3.1. 
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